library(foreign)
library(sandwich)
library(stargazer)
setwd("~/Replication Code/")

##---Loading in data ----
anes9297 <- read.dta("./data/anes9297_subset.dta")

ccap08 <- read.dta("./data/ccap08_subset.dta")

ccap12 <- read.dta("./data/ccap12_subset.dta")

vsg.dat <- read.dta("./data/vsg_subset.dta")

ccap16 <- read.dta("./data/ccap16_subset.dta")

#----------------------------------------------------------------------------------------------------#
# ANES MODELS  -- 92-94
#----------------------------------------------------------------------------------------------------#
dat1 <- na.omit(anes9297[c("white_int92", "white_int94",
                            "pid92_7scaled", "pid94_7scaled",
                            "rr92_scaled", "rr94_scaled", "V940005")])

### PID on RR
m1.std <- lm(scale(rr94_scaled) ~  scale(pid92_7scaled)  + scale(rr92_scaled), 
             data = dat1, weights = V940005, 
             subset = white_int92 == white_int94)
summary(m1.std)
rse.m1.std <- sqrt(diag(vcovHC(m1.std, type = "HC1")))

### RR on PID
m10.std <- lm(scale(pid94_7scaled) ~  scale(pid92_7scaled) + scale(rr92_scaled),
              data = dat1, weights = V940005, 
              subset = white_int92 == white_int94)
summary(m10.std)
rse.m10.std <- sqrt(diag(vcovHC(m10.std, type = "HC1")))

#----------------------------------------------------------------------------------------------------#
# CCAP 2008
#----------------------------------------------------------------------------------------------------#
dat2 <- na.omit(ccap08[c("pid7r_m_sc", "rr_m_scaled", "rr_o_scaled", 
                          "pid7r_o_sc", "weight")]) 

### PID on RR
m2.std <- lm(scale(rr_o_scaled) ~ scale(pid7r_m_sc)  + scale(rr_m_scaled), 
             data = dat2, weights = weight)
summary(m2.std)
rse.m2.std <- sqrt(diag(vcovHC(m2.std, type = "HC1")))


### RR on PID
m20.std <- lm(scale(pid7r_o_sc) ~ scale(pid7r_m_sc)  + scale(rr_m_scaled), 
              data = dat2, weights = weight)
summary(m20.std)
rse.m20.std <- sqrt(diag(vcovHC(m20.std, type = "HC1")))

#----------------------------------------------------------------------------------------------------#
# CCAP 2012: March
#----------------------------------------------------------------------------------------------------#
ccap.p2 <- subset(ccap12, p2 == 1, 
                  select=c("rr_p1_scaled", 
                           "rr_p2_scaled", "pp_pid7r_sc", "pid7_sc", "weight"))
dat3 <- na.omit(ccap.p2)
# PID on RR

m3.std <- lm(scale(rr_p2_scaled) ~  scale(pp_pid7r_sc)  + scale(rr_p1_scaled), 
             data = dat3, weights = weight)
summary(m3.std)
rse.m3.std <- sqrt(diag(vcovHC(m3.std, type = "HC1")))


# RR on PID
m30.std <- lm(scale(pid7_sc) ~  scale(pp_pid7r_sc)  + scale(rr_p1_scaled), 
              data = dat3, weights = weight)
summary(m30.std)
rse.m30.std <- sqrt(diag(vcovHC(m30.std, type = "HC1")))

#----------------------------------------------------------------------------------------------------#
# CCAP 2012: August
#----------------------------------------------------------------------------------------------------#
ccap.p3 <- subset(ccap12, p3 == 1, 
                  select=c("rr_p1_scaled", "rr_p3_scaled", "pp_pid7r_sc", "pid7_sc", "weight"))
dat4 <- na.omit(ccap.p3)
# PID on RR
m4.std <- lm(scale(rr_p3_scaled) ~  scale(pp_pid7r_sc)  + scale(rr_p1_scaled), 
             data = dat4, weights = weight)
summary(m4.std)
rse.m4.std <- sqrt(diag(vcovHC(m4.std, type = "HC1")))


# RR on PID
m40.std <- lm(scale(pid7_sc) ~  scale(pp_pid7r_sc)  + scale(rr_p1_scaled), 
              data = dat4, weights = weight)
summary(m40.std)
rse.m40.std <- sqrt(diag(vcovHC(m40.std, type = "HC1")))

#----------------------------------------------------------------------------------------------------#
# VSG 2012-2016
#----------------------------------------------------------------------------------------------------#
dat5 <- na.omit(vsg.dat[c("rr_sc_12", "rr_sc_16", "pid7_sc_12", "pid7_sc_16", "weight")])

# PID on RR
m5.std <- lm(scale(rr_sc_16) ~  scale(pid7_sc_12)  + scale(rr_sc_12), 
             data = dat5, weights = weight)
summary(m5.std)
rse.m5.std <- sqrt(diag(vcovHC(m5.std, type = "HC1")))

# RR on PID
m50.std <- lm(scale(pid7_sc_16) ~  scale(pid7_sc_12)  + scale(rr_sc_12), 
              data = dat5, weights = weight)
summary(m50.std)
rse.m50.std <- sqrt(diag(vcovHC(m50.std, type = "HC1")))

#----------------------------------------------------------------------------------------------------#
# CCAP 2016
#----------------------------------------------------------------------------------------------------#
dat6 <- na.omit(ccap16[c("b_rr_sc", "p_rr_sc", "b_pid7_sc", "p_pid7_sc", "weight_post")])

# PID on RR
m6.std <- lm(scale(p_rr_sc) ~  scale(b_pid7_sc)  + scale(b_rr_sc), 
             data = dat6, weights = weight_post)
summary(m6.std)
rse.m6.std <- sqrt(diag(vcovHC(m6.std, type = "HC1")))

# RR on PID
m60.std <- lm(scale(p_pid7_sc) ~  scale(b_pid7_sc)  + scale(b_rr_sc), 
              data = dat6, weights = weight_post)
summary(m60.std)
rse.m60.std <- sqrt(diag(vcovHC(m60.std, type = "HC1")))

# Table -------
m1.std_p <- m1.std 
m10.std_p <- m10.std 
m2.std_p <- m2.std 
m20.std_p <- m20.std 
m3.std_p <- m3.std 
m30.std_p <- m30.std
m4.std_p <- m4.std 
m40.std_p <- m40.std 
m5.std_p <- m5.std 
m50.std_p <- m50.std 
m6.std_p <- m6.std 
m60.std_p <- m60.std

names(m1.std_p$coefficients) <- names(coef(m1.std))
names(m10.std_p$coefficients) <- names(coef(m1.std))
names(m2.std_p$coefficients) <- names(coef(m1.std))
names(m20.std_p$coefficients) <- names(coef(m1.std))
names(m3.std_p$coefficients) <- names(coef(m1.std))
names(m30.std_p$coefficients) <- names(coef(m1.std))
names(m4.std_p$coefficients) <- names(coef(m1.std))
names(m40.std_p$coefficients) <- names(coef(m1.std))
names(m5.std_p$coefficients) <- names(coef(m1.std))
names(m50.std_p$coefficients) <- names(coef(m1.std))
names(m6.std_p$coefficients) <- names(coef(m1.std))
names(m60.std_p$coefficients) <- names(coef(m1.std))

names(rse.m1.std) <- names(coef(m1.std)) 
names(rse.m10.std) <- names(coef(m1.std)) 
names(rse.m2.std) <- names(coef(m1.std)) 
names(rse.m20.std) <- names(coef(m1.std))
names(rse.m3.std) <- names(coef(m1.std)) 
names(rse.m30.std) <- names(coef(m1.std)) 
names(rse.m4.std) <- names(coef(m1.std)) 
names(rse.m40.std) <- names(coef(m1.std))
names(rse.m5.std) <- names(coef(m1.std)) 
names(rse.m50.std) <- names(coef(m1.std)) 
names(rse.m6.std) <- names(coef(m1.std)) 
names(rse.m60.std) <- names(coef(m1.std))

stargazer(m1.std_p, m10.std_p, m2.std_p, m20.std_p, m3.std_p, m30.std_p, m4.std_p, m40.std_p, m5.std_p, m50.std_p, m6.std_p, m60.std_p, 
          title = "Relationship between Partisanship and Racial Attitudes",
          covariate.labels = c("Partisanship$_{t-1}$",
                               "Racial Resentment$_{t-1}$"),
          dep.var.labels = c("Racial Resentment", "Partisanship", "Racial Resentment", "Partisanship"), 
          se = list(rse.m1.std, rse.m10.std, rse.m2.std, rse.m20.std, 
                    rse.m3.std, rse.m30.std, rse.m4.std, rse.m40.std, 
                    rse.m5.std, rse.m50.std, rse.m6.std, rse.m60.std),
          no.space = T, notes = c("OLS regression results. Robust standard errors in parentheses. All variables standardized."), 
          notes.align = "l", intercept.bottom = T,
          label = c("std_vars"),
          digits = 3, df = F, omit.stat = c("f", "adj.rsq"), align = T, star.char = c("*"), star.cutoffs = c(0.05)
)


#----------------------------------------------------------------------------------------------------#
### Affect Attitude operationalization (Table D.2)
#----------------------------------------------------------------------------------------------------#
## ANES
dat1 <- na.omit(anes9297[c("pid92_7scaled", "pid94_7scaled", "V940005",
                            "white_int92","white_int94",
                            "ftdif_wb92_sc", "ftdif_wb94_sc")])
# PID on Race
m1_1.std <- lm(scale(ftdif_wb94_sc) ~  scale(pid92_7scaled)  + scale(ftdif_wb92_sc), 
           data = dat1, weights = V940005,
           subset = white_int92 == white_int94)
summary(m1_1.std)
rse.m1_1.std <- sqrt(diag(vcovHC(m1_1.std, type = "HC1")))
# Race on PID
m10_1.std <- lm(scale(pid94_7scaled) ~  scale(pid92_7scaled)  + scale(ftdif_wb92_sc), 
            data = dat1, weights = V940005,
            subset = white_int92 == white_int94)
summary(m10_1.std)
rse.m10_1.std <- sqrt(diag(vcovHC(m10_1.std, type = "HC1")))


# VSG 
dat2 <- na.omit(vsg.dat[c("nwa_sc_12", "nwa_sc_16",
                           "pid7_sc_12", "pid7_sc_16", "weight")])

# PID on affect
m2_1.std <- lm(scale(nwa_sc_16) ~ scale(pid7_sc_12) + scale(nwa_sc_12), 
             data = dat2, weights = weight)
summary(m2_1.std)
rse.m2_1.std <- sqrt(diag(vcovHC(m2_1.std, type = "HC1")))

# affect on PID
m20_1.std <- lm(scale(pid7_sc_16) ~ scale(pid7_sc_12) + scale(nwa_sc_12), 
              data = dat2, weights = weight)
summary(m20_1.std)
rse.m20_1.std <- sqrt(diag(vcovHC(m20_1.std, type = "HC1")))


# CCAP 2016
dat3 <- na.omit(ccap16[c("b_pid7_sc", "p_pid7_sc", "weight_post", "b_wb_fav_dif", "p_wb_fav_dif")])


# PID on Net Fav
m3_1.std <- lm(scale(p_wb_fav_dif) ~  scale(b_pid7_sc)  + scale(b_wb_fav_dif), 
             data = dat3, weights = weight_post)
summary(m3_1.std)
rse.m3_1.std <- sqrt(diag(vcovHC(m3_1.std, type = "HC1")))


# Net Fav on PID
m30_1.std <- lm(scale(p_pid7_sc) ~  scale(b_pid7_sc)  + scale(b_wb_fav_dif), 
              data = dat3, weights = weight_post)
summary(m30_1.std)
rse.m30_1.std <- sqrt(diag(vcovHC(m30_1.std, type = "HC1")))

# Table -------
m1_1.std_p <- m1_1.std 
m10_1.std_p <- m10_1.std 
m2_1.std_p <- m2_1.std 
m20_1.std_p <- m20_1.std 
m3_1.std_p <- m3_1.std 
m30_1.std_p <- m30_1.std


names(m1_1.std_p$coefficients) <- names(coef(m1_1.std))
names(m10_1.std_p$coefficients) <- names(coef(m1_1.std))
names(m2_1.std_p$coefficients) <- names(coef(m1_1.std))
names(m20_1.std_p$coefficients) <- names(coef(m1_1.std))
names(m3_1.std_p$coefficients) <- names(coef(m1_1.std))
names(m30_1.std_p$coefficients) <- names(coef(m1_1.std))

rse_1.m1_1.std <- names(coef(m1_1.std)) 
rse_1.m10_1.std <- names(coef(m1_1.std)) 
rse_1.m2_1.std <- names(coef(m1_1.std)) 
rse_1.m20_1.std <- names(coef(m1_1.std))
rse_1.m3_1.std <- names(coef(m1_1.std)) 
rse_1.m30_1.std <- names(coef(m1_1.std)) 

stargazer(m1_1.std_p, m10_1.std_p, m2_1.std_p, m20_1.std_p, m3_1.std_p, m30_1.std_p, 
          title = "Relationship between Partisanship and Affect Differential, Standardized Variables",
          covariate.labels = c("Partisanship$_{t-1}$",
                               "Affect Difference$_{t-1}$"),
          dep.var.labels = c("Affect Difference", "Partisanship", "Affect Difference", "Partisanship"), 
          se = list(rse.m1_1.std, rse.m10_1.std, rse.m2_1.std, rse.m20_1.std, 
                    rse.m3_1.std, rse.m30_1.std),
          no.space = T, notes = c("OLS regression results. Robust standard errors in parentheses. All variables standardized."), 
          notes.align = "l", intercept.bottom = T,
          digits = 3, df = F, omit.stat = c("f", "adj.rsq"), align = T, star.char = c("*"), star.cutoffs = c(0.05)
)

